﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary.NumericalIntegration
{
    /// <summary>
    /// This class provides a bivariate numerical integration by using the double exponential transform method.
    /// </summary>
    [Serializable]
    public class DoubleExponentialTransformIntegrator : AbstractIntegrator
    {
        /// <summary>
        /// Initializes a new instance of the <see cref="DoubleExponentialTransformIntegrator"/> class.
        /// </summary>
        /// <param name="a">The lower value a of the integral.</param>
        /// <param name="b">The upper value b of the integral.</param>
        /// <param name="function">The function of the integral.</param>
        public DoubleExponentialTransformIntegrator(double a, double b, IRealFunction function)
            : base(a, b, function)
        {
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="DoubleExponentialTransformIntegrator"/> class.
        /// </summary>
        /// <param name="a">The lower value a of the integral.</param>
        /// <param name="b">The upper value b of the integral.</param>
        /// <param name="function">The function of the integral.</param>
        public DoubleExponentialTransformIntegrator(double a, double b, SimplePolynomial function)
            : base(a, b, function)
        {
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="DoubleExponentialTransformIntegrator"/> class.
        /// </summary>
        /// <param name="a">The lower value a of the integral.</param>
        /// <param name="b">The upper value b of the integral.</param>
        /// <param name="function">The function of the integral.</param>
        public DoubleExponentialTransformIntegrator(double a, double b, Polynomial function)
            : base(a, b, function)
        {
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="DoubleExponentialTransformIntegrator"/> class.
        /// </summary>
        /// <param name="a">The lower value a of the integral.</param>
        /// <param name="b">The upper value b of the integral.</param>
        /// <param name="function">The function of the integral.</param>
        public DoubleExponentialTransformIntegrator(double a, double b, SimpleRationalPolynomial function)
            : base(a, b, function)
        {
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="DoubleExponentialTransformIntegrator"/> class.
        /// </summary>
        /// <param name="a">The lower value a of the integral.</param>
        /// <param name="b">The upper value b of the integral.</param>
        /// <param name="function">The function of the integral.</param>
        public DoubleExponentialTransformIntegrator(double a, double b, HardRealFunction function)
            : base(a, b, function)
        {
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="DoubleExponentialTransformIntegrator"/> class.
        /// </summary>
        /// <param name="integral">The integral for the numerical integration.</param>
        public DoubleExponentialTransformIntegrator(FinitelyDoubleIntegral integral)
            : base(integral.A, integral.B, integral.Function)
        {
        }

        /// <summary>
        /// Integrates the integral in a numerical way.
        /// </summary>
        /// <param name="iterations">The maximum numbers of iterations.</param>
        /// <returns>The area of the specified integral.</returns>
        public double Integrate(int iterations)
        {
            double a = this.A;
            double b = this.B;
            bool boundSwap = false;

            if (b < a)
            {
                double tempuri = this.A;

                a = b;
                b = tempuri;
                boundSwap = true;
            }

            double result = InternalIntegrate(this.Function, iterations, a, b);

            if (boundSwap)
            {
                return -result;
            }

            return result;
        }

        private static double InternalIntegrate(IRealFunction f, int n, double a, double b)
        {
            int k, l;
            double z, exz, hcos, hsin, s, dxdz, x1, x2, w;

            n = n / 2;

            double h = 5.0 / n;
            double ss = 0.5 * f.SolveAt((a + b) / 2);
            double ss1 = 0.0;
            double ss2 = 0.0;
            const double pi = Math.PI / 4;

            for (k = -n; k < 0; k++)
            {
                z = h * (double) k;
                exz = Math.Exp(z);
                hcos = exz + 1.0 / exz;
                hsin = exz - 1.0 / exz;
                s = Math.Exp(pi * hsin);
                w = s + 1.0 / s;
                dxdz = hcos / (w * w);
                x1 = (b * s + a / s) / w;
                x2 = (a * s + b / s) / w;

                if (x1 != a && x1 != b)
                {
                    ss1 += dxdz * f.SolveAt(x1);
                }

                if (x2 != a && x2 != b)
                {
                    ss2 += dxdz * f.SolveAt(x2);
                }
            }

            for (l = 0; l < 3; l++)
            {
                for (k = -n; k < 0; k++)
                {
                    z = h * ((double) k + 0.5);
                    exz = Math.Exp(z);
                    hcos = exz + 1.0 / exz;
                    hsin = exz - 1.0 / exz;
                    s = Math.Exp(pi * hsin);
                    w = s + 1.0 / s;
                    dxdz = hcos / (w * w);
                    x1 = (b * s + a / s) / w;
                    x2 = (a * s + b / s) / w;

                    if (x1 != a && x1 != b)
                    {
                        ss1 += dxdz * f.SolveAt(x1);
                    }

                    if (x2 != a && x2 != b)
                    {
                        ss2 += dxdz * f.SolveAt(x2);
                    }
                }

                n *= 2;
                h /= 2;
            }

            return 2 * (b - a) * pi * h * (ss + ss1 + ss2);
        }
    }
}